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Abstract 

We study the cosmological evolution of domain wall networks in two and three 
spatial dimensions in the radiation and matter eras using a large number of high- 
resolution field theory simulations with a large dynamical range. We investigate 
the dependence of the uncertainty in key parameters characterising the evolution 
of the network on the size, dynamical range and number of spatial dimensions 
of the simulations and show that the analytic prediction compares well with the 
simulation results. We find that there is ample evidence from the simulations of a 
slow approach of domain wall networks towards a linear scaling solution. However, 
while at early times the uncertainty in the value of the scaling exponent is small 
enough for deviations from the scaling solution to be measured, at late times the 
error bars are much larger and no strong deviations from the scaling solution are 
found. 
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1 Introduction 


Topological defects are generic in nature and may be formed whenever a phase 
transition occurs. In an expanding universe cooling down from a very hot ini- 
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tial state it is to be expected that topological defects may provide a unique 
window onto the physics of the early universe offering perhaps the best hope of 
a clear observable link between cosmology and particle physics [1,2]. Most cos¬ 
mological studies of cosmic defects have focused on cosmic strings due to their 
interesting properties and strong motivation from fundamental physics (see 
for example [3,4,5] and references therein). Although standard cosmic strings 
are now ruled out as the sole contribution to the large scale structure of the 
universe [6,7] they may still be the dominant source of perturbations on small 
cosmological scales and may give rise to a number of interesting cosmological 
consequences [8,9,10]. However, at present there are only two observations for 
which cosmic strings seem to offer the most natural explanation [11,12]. 

Domain wall scenarios have attracted less attention since heavy domain walls 
in a linear scaling regime rapidly dominate the energy density of the universe. 
Moreover, domain walls which are light enough to satisfy current CMB con¬ 
straints have a negligible direct contribution to structure formation. However, 
in this case a number of interesting consequences are possible such as a con¬ 
tribution to the dark energy [13,14] (if domain walls are frozen in comoving 
coordinates) and a small bnt measurable contribution to the CMB anisotropies 
at large angular scales [15]. Domain walls may also separate regions in the nni- 
verse with different values of the cosmological parameters and/or fundamental 
constants of natnre [16,17]. 

In this paper we perform a quantitative stndy of the cosmological evolution 
of domain wall networks. We investigate the dependence of the uncertainty 
in key parameters characterising the evolution of the network on the size, 
dynamical range and nnmber of spatial dimensions of the simulations using 
a simple analytic model. We then compare onr analytic predictions with the 
results of a large set of high-resolution simnlations of domain walls in two 
and three spatial dimensions [18], using the standard Press-Ryden-Spergel 
(PRS) algorithm [19] (see also [20,21,22,23]), and discuss the evidence from 
the simnlations of a slow approach towards a linear scaling regime. Previous 
studies of domain wall network evolution [19,20,21,22,23] having a smaller 
number of simulations with smaller size and dynamical range than the present 
one have found some hints for deviations from a scale-invariant evolution. It 
is therefore crucial to investigate if these are only transitory or if there is a 
more fnndamental reason for snch deviations. 

The present article is a follow-np of [18]. There, we concentrated on the overall 
(global) dynamical features of the simulations. On the other hand, having a 
large dynamic range means that a more localised analysis is also possible, and 
in particular local exponents can be calculated with relatively small errors. In 
the present paper we explore this possiblity, and also make use of the large 
number of simulations to discuss some analytic ways to estimate statistical 
errors. 
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2 Domain Wall Network Evolution 


We study the evolution of a domain wall network in a flat homogeneous and 
isotropic Friedmann-Robertson-Walker (FRW) universe. We consider a scalar 
held 0 with the Lagrangian density 

^ - ^( 0 ) , ( 1 ) 


and we will take V{<p) to be the generic 0^ potential with two degenerate 
minima given by 


1/(0) = Ro 



( 2 ) 


which obviously admits domain wall solutions. Following the procedure de¬ 
scribed in ref. [19] we modihed the equations of motion in such a way that 
the co-moving thickness of the domain walls is hxed in co-moving coordinates 
allowing us to resolve the domain walls throughout the full dynamical range of 
the simulations. With this modihcation implemented the equations of motion 
for the held 0 become: 


(9^0 f dlna\ d4> ^ 2 , b^V 

^ ^ “ [dlnri 


( 3 ) 


where a is the scale factor, rj is the conformal time and a and (3 are constants. 
We take 0 = 0 in order to have constant co-moving thickness and a = 3 
to ensure that the momentum conservation law of the wall evolution in an 
expanding universe is maintained [19]. Equation (3) is then integrated using 
a standard hnite-diherence scheme. 

We have verihed that the PRS alghoritm gives the correct results in some 
special cases such as the dynamics of a plane wall or the collapse of a spherical 
or cilindrical domain wall. We have also verihed that it appears to have a small 
impact on the large-scale dynamics of domain wall networks and does not seem 
to ahect the quantities we want to measure for the purpose of testing scaling 
properties provided a minimum acceptable tickness is used. However, it is only 
possible to test the performance of the PRS alghoritm over a narrow window 
since the ‘true’ equation of motions for the domain walls rapidly make the 
wall thickness smaller than the grid size. 

In addition to these simple tests with domain walls, the PRS algorithm has 
been much more extensively used and tested in the case of cosmic strings 
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(see for example ref. [24]). From those one can infer that with the use of the 
PRS algorithm some quantitative features of the networks will indeed differ 
(for example, the distribution of small-scale features on the networks), but 
the broad features will be largely unchanged. An example of the latter is the 
existence of an attractor scaling solution, and how fast it is reached starting 
from some given initial conhguration. Since this is the issue we are studying 
here, we believe that the PRS algorithm is adequate for our purposes. 

The ratio between the kinetic and potential energy of the domain walls is 
approximately given by 



where A is the total co-moving area of the domain walls determined using the 
algorithm described in ref. [19,22,18] and we are measuring length in units 
of the grid spacing Ax (so that Ax = 1). This quantity is related to the 
root-mean squared velocity of the domain walls which should be conserved 
in a linear scaling regime. We assume the initial value of 0 to be a random 
variable between —00 and 4-00 and the initial value of d(j)/dri to be zero. See 
[22,18,24] for further discussion of these and other issues. 


3 Analytic Modelling 


We consider a simple model for the evolution of the uncertainty in key param¬ 
eters characterising the evolution of a domain wall network which we then test 
against domain wall network simulations. Let us dehne the comoving correla¬ 
tion length of the network ^ = V/A. Consider a cubic grid in Nd dimensions 
where each cube of comoving volume has Nd faces of its own. Note 

that the number of faces of a cube in dimensions is 2N£) but each one of 
them belongs to two adjacent cubes so that in practise each cube has Nj^ faces 
of its own. We also assume that there is on average one domain wall of comov¬ 
ing area pgj; cube of comoving volume occupying one of the faces of 

the cube. Hence, the probability that a face of such a cube is occupied by a 
domain wall is p = l/N^ so that the variance in the number of domain walls 
per cube is (t| = Ncpq where q = 1 — p. Hence If we now 

consider a cube of comoving volume V = it will contain N = V/V^ cubes 
of comoving volume V^. The standard deviation, ax, of X = n/n (where n 
denotes the number of domain walls of comoving area n its average 
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value) on a given volume, V, is proportional to y/JV/N and consequently 


<^x = o-^N 




-Nol2 



( 5 ) 


This can also be calculated for any other variable, Y, characterising the net¬ 
work if Y is proportional to X. In this case 


ay = (JxY. 


( 6 ) 


Here the random variable Y may represent the total area A, or the ratio 
between the kinetic and potential energy parametrised by F. Although in the 
hrst case it is a fair assumption to draw a direct relation between the number 
of domain walls and the total area, in the second case we expect the relation 
to be less direct due to the dispersion in the domain wall velocities and to the 
existence of an important fraction of the kinetic energy which is not associated 
with the domain walls themselves but with their decay products. 


The scaling exponent may be calculated from the value of R = rjA/V at two 
different values of the conformal time rji and 772 as 


_ In {R1/R2) In (^1/^2) Ai - A2 

- - rs_/ - —|— - 

In (771/772) In (771/772) In (771/772) 


( 7 ) 


so that 


In ( 771 / 772 ) 


( 8 ) 


assuming that Ai and A 2 are small (compared with .Ri and R 2 respectively) 
and uncorrelated. Here A = {R — R)/R where R{ri) denotes the average value 
of R{f]) over a given number. Ns, of simulations. We also may also estimate 
ay for several key parameters, Y = R,F,X, describing the evolution of the 
network directly from the simulations as 


Sy = 


Ns-I 


Ns 

EXi 


Yf 


i=l 


( 9 ) 


Note that the average value of Y for a given sample of simulations, Y, is the 
best estimator of the real average of Y with and error of ay/y/N 5 . Also sy is 
the best estimator of ay and sy ~ ay if Ns is large enough. 
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Fig. 1. The evolution of R = rjA/V and F for one hundred 4096^ 2D matter 
era simulations. We also plot the average evolution of these parameters and the 
expected 1-sigma interval around the mean (using the analytic estimates for an 
and ap described in Sec. III). 




Fig. 2. The same as Fig. 1 but now for one hundred 4096^ 2D radiation era simu¬ 
lations. 

4 Results and discussion 


A total of several thousand matter and radiation era simulations in two and 
three spatial dimensions were run for various box sizes and dynamical ranges. 
Here, we compare some of these numerical results with analytic expectations 
and discuss the main results. In Fig. 1 we plot the evolution of R = r]A/V 
and F for one hundred 2D matter era simulations. We also plot the average 
evolution of these parameters and the expected 1-sigma interval around the 
mean (using the analytic estimates for ap and ap described in the previous 
section). Fig. 2 shows a similar plot but now for one hundred 4096^ 2D radi¬ 
ation era simulations. Figs. 3 and 4 display analogous results for one hundred 
256^ 3D matter and radiation era simulations respectivelly. 

To zeroth order the results do not seem to be very sensitive on the number of 
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Fig. 3. The same as Fig. 1 but now for one hundred 256^ 2>D matter era simulations. 



Fig. 4. The same as Fig. 1 but now for one hundred 256^ 2>D radiation era simula¬ 
tions. 

spatial dimensions, No- We clearly see that althongli the network is not scaling 
dnring most of the simnlation it seems to approach a scaling solntion very 
slowly at late times. Also, the main difference between matter and radiation 
era simnlations lies in the larger valne of the typical velocity of the domain 
walls in the latter case althongh we also hnd a considerable difference in the 
evolntion of the matter and radiation era rnns at early times when the initial 
conditions are still important. 

As expected we see that in the case of the domain wall area the analytic 
estimate is a good approximation {aR is only slightly overestimated - see Figs. 
7 and 8) while in the case of the parameter F the nncertainty is nnderestimated 
dne to not taking into acconnt fact the nncertainty in the velocity of the 
domain walls as well as the part of the kinetic energy which is not associated 
with the domain walls themselves bnt with the particle radiation emitted by 
them. 

Onr analytic estimates for the variance of the scaling exponent. A, were based 
on the assnmption that the statistical properties of the domain wall network 
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Fig. 5. Comparison between numerical estimates of the standard deviation of the 
scaling parameter in an ensemble of Ns simulations, sx/V^S^ with sx computed 
either directly from the simulation (inner bars) or calculating sa from the simu¬ 
lations and using eqn. ( 8 ) to compute sx assuming no correlation (outer bars) for 
fli = f]* = and various values of r ]2 = rj for one hundred 4096^ 2D and one 

hundred 256^ 3D matter era simulations (left and right plots respectivelly). 
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Fig. 6 . The same a Fig. 5 but now for one hundred 4096^ 2D and one hundred 
256^ 3D radiation era simulations (left and right plots respectivelly). As expected 
when we increase r ]2 — r]i our assumption of uncorrelation becomes increasingly more 
justified. The error bars are always placed in the middle of the interval [rii^'q 2 \- 

at two different times r/i and 772 were nncorrelated. Here we test for the validity 
of that assumption for different sizes of the conformal time interval 772 — h- 
In Fig. 5 we compare the numerical estimates of the standard deviation of 
the scaling parameter in an ensemble of Ns simulations, sx/V^s, with sx 
computed either directly from the simulation (inner bars) or calculating sa 
from the simulations and using eqn. (8) to compute Sx assuming no correlation 
(outer bars) for 771 = 77 * = and various values of 772 = 77 for one 

hundred 4096^ 2D and one hundred 256^ 3D matter era simulations (left and 
right plots respectivelly). We clearly see that when we increase 772 — 771 our 
assumption of no correlation becomes increasingly more justihed and that 
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Fig. 7. Evolution of the scaling exponent A as a function of the conformal time r] 
for one hundred 4096^ 2D and one hundred 256^ matter era simulations (left 
and right plots respectivelly). The binnings have constant dynamical range (with 
— 2) and the error bars are the standard deviation of the scaling parameter 
in an ensemble of Ng simulations, sx/VNg, calculated either using our analytic 
estimate for a\ (outer bars) or by calculating sa from the simulations and using eqn. 
(8) to compute s\ assuming no correlation (inner bars). All the error bars (except 
the last one in each plot) were artificially enlarged by a factor of where 

rjf is the conformal time corresponding to the last error bar. 



log(ii) 


Fig. 8. The same a Fig. 7 but now for one hundred 4096^ 2D and one hundred 256^ 
3D radiation era simulations (left and right plots respectivelly). 

there is no significant dependence of the correlation time on the nnmber of 
spatial dimensions No- Fig. 6 displays analogous results but now for radiation 
era simulations. However, no signihcant differences from the matter era runs 
are found. 


In Fig. 7 we plot the evolution of the scaling exponent A as a function of the 
conformal time r] for one hundred 4096^ 2D and one hundred 256^ SD matter 
era simulations (left and right plots respectivelly). Here, we chose binnings of 
constant dynamical range (with r] 2 /r]i = 2) and the error bars are the standard 
deviation of the scaling parameter in an ensemble of Ng simulations, s\/\/Ng, 
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calculated either using our analytic estimate for a\ (outer bars) or by calcu¬ 
lating sa from the simulations and using eqn (8) to compute s\ assuming no 
correlation (inner bars). Note that the error bars in Fig. 7 were artificially 
enlarged by a factor of of where rjf is the conformal time cor¬ 

responding to the last error bar on the right so that at early times the real 
error bars are much smaller than at late times. We clearly see that our simple 
model successfully predicts the evolution of the uncertainties in the value of 
the scaling exponent, A, although the performance of the method may slightly 
depend on the number of spatial dimensions, N^,. We hnd that, as expected, 
at early times there are strong deviations from the linear scaling solution but 
the networks seem to slowly approach a linear scaling solution at late times. 
Fig. 8 displays analogous results but now for radiation era simulations but 
no signihcant differences (as far as the approach to scaling is concerned) were 
found. 

Our results are consistent with those of previous studies [19,20,21,22,23] in 
hnding a slow approach to linear scaling as well as strong deviations from a 
linear scaling solution at early times. However, we do not hnd any evidence for 
strong deviations from a linear scaling solution at late times. We note that this 
is not in disagreement with the results obtained by other authors but it is a 
consequence of the larger number, size and dynamical range of the simulations 
analysed in our paper. 

We have also performed a small number of 128^ 4D simulations which in a sim¬ 
ple phenomenologogical way be of interest to brane world scenarios [18,25,26]. 
The results obtained are consistent with those described in this paper for the 
2D and 3iA simulations. 
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